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Abstract 



Shell corrections of finite, spherical, one-body potentials are analyzed us- 
ing a smoothing procedure which properly accounts for the contribution from 
the particle continuum, i.e., unbound states. Since the plateau condition for 
the smoothed single-particle energy seldom holds, a new recipe is suggested for 
the definition of the shell correction. The generalized Strutinsky smoothing 
procedure is compared with the results of the semi-classical Wigner-Kirkwood 
expansion. A good agreement has been found for weakly bound nuclei in 
the vicinity of the proton drip line. However, some deviations remain for 
extremely neutron-rich systems due to the pathological behavior of the semi- 
classical level density around the particle threshold. 
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I. INTRODUCTION 



Positive-energy eigenstates of the average single-particle potential are very important for 
the description of nuclei close to the particle drip lines where the Fermi level approaches 
zero (see Ref. |l| and references quoted therein), and in analysis of highly excited nuclear 
modes such as giant resonances [@,[|. With the advent of radioactive nuclear beams, of 
particular interest are masses of weakly bound nuclei with extreme N/Z ratios. For these 
nuclei, important for both nuclear structure and astrophysics, special care should be taken 
when dealing with the particle continuum, which strongly influences many nuclear properties, 
including global ones (e.g., masses, radii, shapes) as well as nuclear dynamics (i.e., excitation 
modes). 

In a previous paper [|J, a macroscopic- microscopic method was applied to nuclei far 
from the beta stability line. It has been demonstrated that the systematic error in binding 
energies, due to the particle continuum, can be as large as several MeV at the neutron drip 
line; hence it can seriously affect theoretical mass predictions for nuclei far from stability. 
This error depends weakly on deformation, thus suggesting a possibility of re-normalizing 
potential energy surfaces at the spherical shape. 

In this paper, the effect of the particle continuum on the shell correction (the quantal 
contribution to the total energy in the macroscopic-microscopic approach) is investigated by 
solving the Schrodinger equation in the complex energy plane. The new procedure allows 
for the direct treatment of both narrow resonances and the smooth continuum background 
when calculating the single-particle level density. 

The paper is organized as follows. Section || contains a brief review of the shell-correction 
method in terms of the single-particle level density. The semi-classical approach is discussed 



in Sec. III. Section IV describes the modified Strutinsky re- normalization procedure which 



takes care of the continuum effects. The results of the calculations are contained in Sec. [V]. 



Finally, conclusions are given in Sec. [V| The threshold behavior of the semi-classical level 
density is discussed in Appendix |A]. 

II. SHELL CORRECTION AND SINGLE-PARTICLE LEVEL DENSITY 

In the standard macroscopic- microscopic approach ||-§|, the shell correction 

5E = E s . p . - E s . p . (1) 
is the difference between the total single-particle energy -E^.p., 

E s .p. = £ i ' ( 2 ) 



l — occ 



and the smooth single-particle energy E s . p .. The shell correction represents the fluctuating 
part of the binding energy resulting from the single-particle shell structure. 

For the sake of simplicity, we shall restrict our discussion to spherically symmetric nuclei 
and assume that the single-nucleon energy spectrum is that of a one-body Hamiltonian 
H = T+V with a finite, local, and spherically symmetric potential V(r). Since the spectrum 
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of H contains both bound < 0) and unbound (e > 0) eigenvalues, the single-particle level 
density is the sum of the discrete and continuum contributions |)-13] 



g{e)= J2(23i + l)5(e-e i )+g c (e). (3) 

i 

The continuum part, g c (s), is defined in terms of the scattering phase shifts Sij(e) 

M = IX W + 1) «M. (4) 

1,3 

In the shell- correct ion method |5]|| the smooth energy E s >p> is calculated by employing 
the smooth level density g(e) obtained from g(e) by folding with a smoothing function f(x): 

1 r+°° f e' — e\ 

m = ~ de' g(e')fl £ -^). (5) 



7 j-oo y 7 

The folding function f(x) can be written as a product of a weighting function and a 
curvature correction polynomial of the order p=2M 0. The smoothing procedure should 
be unambiguous, i.e., the averaging should leave the smooth part of the level density un- 
touched (the so-called self-consistency condition for the Strutinsky smoothing). This defines 
a curvature correction polynomial for any specific choice of weighting function. In this study, 
a Gaussian weighting function, ^exp(— x 2 ) has been used. The corresponding curvature 

function is an associated Laguerre polynomial L^J^ix 2 ). This choice guarantees the self- 
consistency condition for g(e) if the smooth level density behaves as a polynomial of degree 
2M + 1 or lower in e. 

The smoothed level density (M) defines both the smooth single-particle energy 



E a . p . = / eg(e)de, (6) 

J — oo 

and the smoothed Fermi level A. The latter is obtained from the particle number equation: 

N = f X g{e)de, (7) 



where N is the total number of particles (i.e. protons or neutrons). The smoothing range 
7 should be greater than the average energy distance between neighboring major shells, 
hn^Al/A 1 / 3 MeV 0. 



Since the result of the smoothing procedure should not depend on the actual form of the 
smoothing function, in particular on the smoothing range 7 and the order p of curvature 
correction, the smooth energy should satisfy the so-called plateau condition: 

For infinite potentials such as a harmonic oscillator or a deformed Nilsson potential, one can 
always find an interval of the smoothing parameters 7 and p in which the smooth energy, 
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hence the shell correction, is practically independent of the values 7 and p [0,0. For 
finite-depth potentials, additional complications arise due to the presence of the continuum 
contribution, Eq. (f|) (see discussion in Refs. and references quoted therein). In most 
calculations applying the shell-correction approach, the continuum is treated approximately 
by using the quasi-bound states, i.e., the states resulting from the diagonalization of a finite 
potential in a large harmonic oscillator basis [^]. However, for light and weakly bound nuclei, 
the plateau condition (BJ) usually does not hold M. 



III. SEMI-CLASSICAL TREATMENT OF SHELL CORRECTION 



A possible alternative to Strutinsky's smoothing procedure is the semi-classical averaging 
based on the Wigner-Kirkwood expansion |TB|-F2TJ . In Ref . |T7J the equivalence between the 
Strutinsky approach and the Wigner-Kirkwood (WK) expansion has been demonstrated. It 
is important to note that this proof assumed that in the Strutinsky approach the plateau 
condition could be fulfilled. 

In the WK expansion, the diagonal part of the the Bloch density [pl, 22]| is 



h 2 p 2 

12M 



W(r) - ^(W(r)) 2 



+ 



(9) 



The spatial density is obtained by using the inverse Laplace transform 

'C(r,0)\ 



p(r,e) 



zr 1 

P — >s 



ft 



and the particle number integral is given by 

N(e) = J p(r,e) d 3 r. 



(10) 



(11) 



Keeping only the two leading terms in the curly bracket in Eq. (^), the WK particle 
number equation can be expressed explicitly in terms of the single-particle potential 



NJe) 



4 /2MV 



r sc (e) 



r 2 \{s-V)* 



h 2 



V 2 V 



32M {e-V)\ 



dr. 



(12) 



The integral in Eq. ([12] ) is cut off at the classical turning point, r sc (e), defined by the 
relation V(r sc ) = e. (For the inclusion of the spin-orbit term see Ref. ||20|| .) The semi- 
classical value of the Fermi energy, A sc , can be determined from the particle number equation 
N SC (X SC ) = N. 

The semi-classical level density is defined as 



dN sc (e) 
de 



(13) 



Here, it is worth reminding that the semi-classical level density is defined only for e > Vq, 
where Vq denotes the bottom of the potential well. That is, g sc {^) = ^ ii e < Vq. An explicit 



expression for g sc (e) in terms of a WK expansion can be found in, e.g., Ref. |23 
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It is to be noted that Eqs. (P-|T3|) are valid for any smooth potential regardless whether 
it is infinite or not. Often, the semi-classical level density is defined in terms of the partition 
function Z{(3) = J C(r,f3) d 3 r. However, there is a difficulty: Z(/3) diverges for finite 
potentials [ffijfl . In practical calculations it is possible to overcome this serious problem by 
an appropriate modification of the potential at large disstances. Formally, this procedure is 
equivalent to placing an infinite potential well (a box with soft walls) at very large distances 
from the classical region 



Since the particle-number integral (pT2[ ) depends only on the potential in the classical 
region, where V <e<VB, the quantity N sc (e) is well defined for energies which are not too 
close to the top of the potential well Vq where the semi-classical expansion breaks down (see 
Appendix |A]). Hence the total number of particles can be calculated without the explicit 
use of the partition function; there is no need at all to put the system into a box in the case 
of a finite potential. 

The total energy of the system of non-interacting particles is 

E = f X eg{e)de = XN - N{e)de. (14) 



By means of Eqs. (|TTD and flltf) , E can be written as 

-i fC(r,(3Y 



E = XN-J C£ x j d 3 r. (15) 

Here it is assumed that the order of the integration with respect to e (and r) and the inverse 
Laplace transform can be interchanged. Using Eq. ([15]), the semi-classical smoothed binding 
energy can be written as 



E SC = \ SC N- l£,l x J^j^)d\ (16) 



where the notation C sc (r, f3) means that only the terms which are not a higher order than h 
are kept in the WK expansion of C(r,/3). As it was discussed in ||20|| , for the determination 
of A sc it is enough to keep the terms of order h~ in Eq. (P) . The explicit expression for the 



smoothed energy is quite lengthy and can be found in Ref. [20 . 



IV. GENERALIZED SHELL-CORRECTION METHOD 

The impact of the particle continuum on shell corrections has been investigated numer- 
ically in Ref. [O] for neutrons in 208 Pb and 298 114 by explicit calculation of the continuum 



part of the level density, Eq. (|J). They have shown that, by taking into account contri- 
butions from the neutron continuum up to ~100 MeV in 208 Pb, the plateau condition @ 
could be met (see Ref. @ for an updated discussion of the continuum contribution in 208 Pb). 
Based on this early exercise, it was generally assumed that the plateau condition could be 
fulfilled for finite potentials provided that the continuum part was included. A systematic 
study of this assumption is given in Sec. |V B. 



The lack of systematic studies using the continuum level density is due to the fact that 
these calculation are quite cumbersome. Except for some special cases, the solution of the 
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radial differential equation can be done only numerically. The calculation of the phase 
shifts and the continuum level density along the real energy axis should be carried out with 
great care. In general, the use of an extremely fine energy step is required in order to 
collect contributions from narrow resonances. Recently, this difficulty has been overcome by 
applying a new method which uses Gamow states [0. A Gamow resonance is a generalized 
eigenstate of the radial Schrodinger equation corresponding to a complex energy eigenvalue 
Wi — Si — iVi (for bound states 1^=0) . The wave function of a Gamow resonance is regular 
at r=0 and has purely outgoing asymptotics with a discrete complex wave number. 

In the new method (see Ref. || for details ), the smoothed level density g(e) can be 
written with the help of the Cauchy theorem as a sum over bound and resonant states, and 
an integral along a contour L in the complex energy plane: 

~9{e) = E / (^p) + J L dw g c {w)f {^\ . (17) 

In Eq. (|l7r), the summation runs over all the bound states and those resonances which are 
above the contour L and below the real energy axis. 

Thanks to the Cauchy theorem, the level density as given by Eq. (|T7|) has only real part 
(the imaginary parts of the two terms in the r.h.s. of Eq. (|TT| ) should cancel out exactly). 
Furtermore, g(e) should be independent of the shape of the contour L. However, since the 
calculations are carried out numerically, the exact cancellation of the imaginary parts is 
always slightly violated and, in addition, the sum of the real parts slightly depends on the 
shape of the contour. For example, if the contour L is extended to the area of the complex 
energy plane with large imaginary energies (Im(w)<-5 MeV) then broad resonances, i.e., 
those with large revalues should be included. As a result, both terms in the r.h.s. of 
Eq. ([T7D would acquire large imaginary parts which would not cancel completely due to 
the numerical errors. Therefore, the best strategy is to choose the contour in such a way 
that it would include only narrow resonances. With this choice the final result is practically 
independent of the shape of the contour, and the imaginary part of g(e) is negligible (it is 
the order of 10~ 4 or less). Three examples of contours L are shown in Fig. [I]. 

With a reasonable choice for the contour L, the Gamow resonances give the major 
contribution from the continuum; the contour integral gives the remaining (small) part. 
From the smoothed level density ([17]) one can determine A and E B p _ using Eqs. (0) and ([|), 
respectively. 

It is worth noting that another, commonly used method of calculating the continuum 
level density is based on the discretization procedure. Here, one assumes that the nucleus 
is placed inside a very large box (cf. the discussion in Sec. |TTTj on the application of the 
semi-classical expansion to finite potentials). Since the properties of the nucleus itself must 
not depend on the box size, one has to re-normalize the level density by subtracting the 
contribution from the free-gas states [13,24-27] . In Refs. [|l|,^7]], the discretization method 



was applied to investigate the accuracy of the semi-classical expression ( |T5j ) for several 
commonly used potentials, and a good agreement was found in all cases. 



V. RESULTS 
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A. Details of calculations 



In the actual calculations, we have used the average Woods-Saxon (WS) potential, which 
contains a central part, a spin-orbit term, and a Coulomb potential for protons. The 
Coulomb potential has been assumed to be that of a charge (Z — l)e distributed with 



the diffused charge density. We employed the set of WS parameters introduced in Ref. p8 
and the charge density form factor was taken in the WS form. (See Ref. @] for details 
pertaining to the single-particle model.) 

The poles of the S'-matrix, i.e. the eigenvalues Wi, have been calculated by solving the 
radial equation numerically using the computer code GAMOW |]29|| . The contour L has been 
chosen to lie far from these poles. This ensures that the energy dependence of phase shifts 
along the path is smooth; hence one can use relatively large energy steps. The phase shifts 
of scattering states along the path L have been calculated by solving the radial equation 
numerically using the unpublished code ZSCAT in which the complex Coulomb routine 
COULCC was used [0. 

As an illustrative example, the distribution of eigenvalues for the stable nucleus 90 Zr 
(neutrons), neutron drip-line nucleus 122 Zr (neutrons), and proton-rich nucleus 180 Pb (pro- 
tons) is shown in Fig. |], together with the contours L used for calculating the level density 
(|T7|). In the case of 90 Zr, there are four poles close to the e = threshold. They are: 
p 1/2 (™=0.2-z0.19 MeV), / 5/2 (w=2.1-i0.3A MeV), i 13/2 (w=3.7-i0.004 MeV), and h 9/2 
(w=3. 9—i0.03 MeV). As seen in Fig. [I], only the Gamow states with the smallest widths, 
i-e., ii3/2 and h 9 / 2 , have been considered explicitly in the level density calculations; a contri- 
bution from the remaining eigenstates has been accounted for by the integral along the path 
L, i.e., by the second term in Eq. (|17|) . For the neutron-rich nucleus 122 Zr, the number of 
near-threshold Gamow states increases. Here, three Gamow resonances: fr/ 2 (w=0.7— z0.02 
MeV), hg/2 (w=3.2-z0.03 MeV), and i 13/2 (u;=4.5-z0.02 MeV) have been used to define the 
contour that includes them in Fig. |l]b. As seen in Fig. [l|(c), the distribution of the proton 
Gamow eigenvalues in 180 Pb is different. Due to the presence of the Coulomb barrier, there 
appear many very narrow resonances even at relatively high energies, i.e. above 10 MeV; 
therefore we have chosen in this case a contour that includes these narrow resonances. 

The energy dependence of neutron phase shifts in 122 Zr and proton phase shifts in 180 Pb 
along the contour L is illustrated in Fig. |2|. Here are shown the real and imaginary parts of 

^cH = -E(2j + i)^H- (is) 

* 1,3 

According to Eq. @, N c can be interpreted as the "continuum particle number" along the 
contour. As expected, the energy dependence of N c along the path is very smooth. It is also 
seen that the imaginary part of N c is small since the contour does not go far from the real 
energy axis. 

The WK calculations in this paper follow those of Ref. [|J except for the treatment of 
the proton average field. Since the combined WS and Coulomb potentials give rise to a 
non-monotonic field, in this study we had to employ a modified treatment of some of the 
terms in the WK expansion according to Ref. |H[ . 
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B. Modified plateau condition 



In order to check the dependence of the smoothed single-particle energy on 7 and p, we 
made systematic calculations of the shell energy by varying these parameters within rea- 
sonable ranges. Figure § shows three typical examples of our analysis of neutron shell cor- 
rections. The plateau condition is satisfied fairly well for the super-heavy nucleus 298 114 184 . 



(This nucleus was previously studied in Ref. |12| .) Here, for each p value, 8E possesses a 



local minimum in 7, and the minimum energy changes little with p. However, in the cases 
of 146 Gd and 90 Zr, it is impossible to assign a definite value to the neutron shell correction. 
The situation for 146 Gd and 90 Zr shown in Fig. |3| should be considered as typical; the plateau 
condition (||) is seldom satisfied. Therefore, we conclude that the proper treatment of the 
continuum level density does not guarantee that the plateau condition is fulfilled. 

In the cases where the plateau condition was approximately satisfied, like in Fig. 0(a), we 
found a strong correlation between the values of 7 and p . In particular, the behavior of g(e) 
as a function of e was found to be very similar for different values of p and 7 corresponding 
to local minima in 5E. Moreover, the dependence of g(e) was found to be almost linear in 
a wide range of e below A. We also checked that for the cases when the plateau condition 
could not be satisfied [see, e.g., Figs. [|(b,c)], the approximate linearity of g{e) was valid. 
It is worthwhile to point out that for the harmonic oscillator potential, the average level 
density behaves as e 2 , while for the finite square well potential, the leading terms behave as 
\fe flnjyrj]. Hence, a local linear behavior of g(e) for a finite WS potential is not unexpected. 
This observation suggests that an alternative recipe for defining shell correction for finite 
potentials, not based on the plateau condition but rather on the behavior of the smooth 
level density, may be possible. 

It is well known that the realistic value of the smoothing parameter has to lie in a certain 
energy interval f23fl . The value of 7 should be large enough to wipe out shell effects in the 



energy range of a typical distance between shells: 

7 > m. (19) 

On the other hand, its upper limit is defined by the number of states considered in the 
calculations, i.e., 

7 < £max - A, (20) 



and by the energy distance between the Fermi level and the bottom of the well [15 



7 < A - V . (21) 

In practice, the optimal value of 7 for a given p is found by the following procedure. First 
we choose the energy interval [£/,£«] which is lying below A and is wider than the average 
shell distance, e.g., e u -Ei>l.5hQ. In this energy interval, we perform the least squares fit 
to the smoothed level density assuming a linear dependence of g(e) on e. The search for 
optimal 7 begins at a small 7 value below HQ where the shell fluctuations are still present, 
and then 7 is gradually increased until the first minimum in SE is found at j p . This 7 P 
corresponds to the smallest value of 7 for a given p that smooths out the shell fluctuations. 
The corresponding shell correction, 5E = SE P , is taken as the optimal shell correction for 
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this p. This procedure is repeated for higher values of p. If variations of 5E P with p are 
small, then the mean value of 5E p represents the shell correction obtained in this modified 
Strutinsky method. The uncertainty of the procedure is given by the r.m.s. error a=a(5E). 

The smoothed level densities g(e), calculated using the above procedure, are displayed 
in Fig. |] for three different values of p=Q, 10, and 14. Since g(e) is practically independent 
of the value of p we amplified the differences by presenting in Fig. f| the ratios of the level 
densities to g(e) belonging to p = 10. 

The shell energy is also practically independent of the value of p. Table [I] displays the 
calculated proton and neutron shell-correction energies and the corresponding r.m.s. errors. 
Since a increases when going to lighter nuclei where the condition does not hold, we 
limited calculations to nuclei heavier than 40 Ca. The calculations were performed for nuclei 
close to the stability valley and for nuclei with extreme N/Z ratios (both neutron-rich and 
proton-rich). It is seen that the r.m.s. error in SE is always less than 250 keV (also for the 
cases such as 90 Zr or 146 Gd where the plateau condition could not be met). Another source 
of theoretical uncertainty lies in the choice of the fitting range [e;,e u ]. In particular, the 
selection of e u plays a role for weakly bound nuclei with very small values of A. In practice, 
in order to guarantee that the fitting region does not overlap with the threshold area, we 
have adopted the value of e u =\ — hQ and ei=e u — 1.5MX In order to estimate the associated 
theoretical error, we have performed a set of calculations for well-bound nuclei assuming 
a larger value of e u , namely e u =X. The average deviation in 5E between the two sets of 
calculations is around 400 keV. We believe that this number represents a fair estimate of the 
uncertainty of our method. 



C. Comparison with the semiclassical method 

As a next step, we performed the detailed comparison of the generalized Strutinsky 
method with the WK expansion. Table | displays the shell correction SE SC = E s p . — E sc 
calculated using the WK method, together with the difference A=5E — 5E SC . In most cases 
A>0. That is, the semi-classical method yields the average single-particle energy E sc which 
is greater than E B _ P _ obtained in the generalized Strutinsky method. The average value of A 
is about 0.45 MeV and the maximal deviation is about 1.8 MeV for neutrons and 0.9 MeV 
for protons. These deviations exhibit large fluctuations with particle number, and they are 
significantly larger than the uncertainty of the generalized Strutinsky smoothing procedure. 
Considering the excellent agreement between the shell energies calculated in the Strutinsky 
and the WK methods obtained previously [p0[| , and the existing proof of the equivalence 
between these two methods [IT]], the large magnitudes of A in Table | seem to be surprising. 



In order to understand this discrepancy, we analyze the behavior of smoothed single- 
particle level densities obtained in both methods. The semi-classical level density has been 
obtained by means of Eq. flT3|), i.e., by calculating the derivative of N sc (e). 

Smoothed level densities g(e) and g sc {£) are compared in Figs. [| and |5] for 146 Gd (neu- 
trons) and 208 Pb (protons), respectively. The behavior of average single-particle densities 
displayed in Fig. |^ shows a generic pattern characteristic of a WS-like potential well [|T3| , |2T 



Namely, g(e) increases monotonically with e reaching the maximum value around the e=0 
threshold for the neutrons (and around Vb, i.e. the top of the Coulomb barrier for the 
protons), and then it smoothly falls down reflecting the increasing contribution from the 
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free-gas states. As discussed in Sec. |VB| , there exists a wide energy region in which g sc (e) 
increases fairly linearly with e. In general, g(e) and g sc (e) are very close except for the 
bottom of the potential well (e~Vo) and close to the top of the potential barrier. 

Considering the low-energy region, there are problems with both methods. The Struti- 
nsky smoothed density is non-zero for e < Vo, i.e., in the classically forbidden region. Here, 
the inequality (|2l|) cannot be met and the averaging method breaks down [15|. Also, there 



are serious questions regarding the applicability of the semi-classical treatment close to the 
s=Vq limit. The fi~ l term in the WK expansion of the semi-classical level density [related 
to the second term in the integral (12)] gives rise to the singularity around the bottom of 
the well. (For a square well potential this singularity behaves as (e — Vo)" 1 / 2 .) As discussed 
in Ref. ||23|| , in the strict treatment of g sc {£), there appears a correction to the level den- 



sity, proportional to the Dirac delta 8{e — Vo). This additional term, usually ignored in 
calculations, partly corrects for a pathological behavior around the bottom of the well by 
introducing a small shift in the Fermi level A sc ||16|| . Table | displays the Fermi energies A 
and A sc . Usually, the Fermi levels are very close; for well-bound nuclei the difference is less 
than 100 keV. Although small, this shift partly contributes to the calculated values of A 
(e.g., for protons in 208 Pb). 

Another major deviation between g(e) and g sc (s) is seen in the energy region close to the 
e=0 threshold in the neutrons, and around the top of the proton Coulomb barrier. Here, 
the reason for this abnormal behavior is our semi-classical approximation. As discussed in 
Appendix 0, the WK neutron level density g sc {£) diverges as ln(— e) / \/—e when e — > 0. For 
protons, the singularity is even more severe: g sc {e) diverges as (Vb — around the top of 
the barrier. 

The pathological behavior of g sc {e) around zero energy is the reason for the largest 
deviations A seen in Table | for the neutron case (e.g., A=1.39MeV for 78 Ni, and it is 
greater than 1 MeV for the Zr isotopes with t4=110,122,124). These nuclei are weakly 
bound (as shown by their low Fermi energies), and the level density around the Fermi level, 
g sc (A), is affected by the threshold effect. 

In order to understand the systematic behavior of A in Table |, neutron shell corrections 
and Fermi energies for the Zr isotopes are shown in Fig. ^ as functions of N. The calculations 
were performed using the single-particle potential corresponding to 90 Zr. The associated 
smoothed level densities are shown in Fig. [7]. It is seen that although the general pattern of 
SE is similar in the generalized shell-correction method and WK approach, A exhibits the 
oscillatory behavior as a function of particle number. The agreement between SE and SE SC 
is very good up to iV~70 (A~-4MeV), but it is spoiled at large neutron numbers where A 
systematically increases. Indeed, above point "C" in Fig. 0, the semi-classical level density 
diverges, and it does not yield a good estimate of the shell correction. 

For the protons, the "dangerous" threshold region of g sc is shifted to higher energies 
due to the presence of the Coulomb barrier (see Fig. [5]) . That explains why the differences 
between the Strutinsky and WK results are smaller for the protons than for the neutrons. 
For instance, for the proton drip-line nucleus 100 Sn the agreement between two methods is 
surprisingly good in spite of the fact that A=0.72 MeV. Since in actual nuclei (both particle- 
bound and in proton emitters) the proton Fermi level is significantly lower than Vb, the 
divergent behavior of proton g sc around the top of the Coulomb barrier has no practical 
importance. 
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VI. CONCLUSIONS 



This paper introduces a new method of calculating the nuclear shell energy. The gen- 
eralized Strutinsky procedure fully takes into account the effect of the particle continuum. 
Although the traditional plateau condition can seldom be met for finite potentials, the pro- 
posed method makes it possible to define shell correction unambiguously. A conservative 
estimate of the uncertainty in SE obtained using the new smoothing procedure is ~400 keV. 
This error can be considered as small. 

In most cases, the results of the generalized Strutinsky procedure are in good agreement 
with those of the semi-classical WK method. Significant deviations have been obtained for 
neutron-rich nuclei for which the neutron Fermi energy is low (A>-4 MeV). This discrepancy 
has been tracked back to the singularity in the WK level density around the top of the 
potential barrier. The density g(e) obtained in the generalized Strutinsky method nicely 
interpolates through the threshold region (see also Refs. I3|.2"7|]). Other advantages of the 



new method are: (i) its applicability to potentials with discontinuous derivatives (e.g., the 
Coulomb potential of a uniform charge distribution and a folded- Yukawa potential) where 
the standard WK expansion cannot be carried out, and (ii) simple generalization to the 
deformed case where the semiclassical expansion becomes awkward ||. 

Finally, let us comment on the difference A between shell corrections obtained in both 
methods. There are several factors which contribute to A. In addition to the threshold 
anomaly mentioned above, other factors are: (i) the shift between the Fermi levels caused 
by the different behavior of the level densities at the bottom of the potential well, (ii) the 
assumption of the local linearity of the smoothed level density, and (iii) the systematic errors 
accumulated during numerical calculations. For protons, our calculations give |A|<900keV 
in all cases considered. Here, the main source of the difference is the deviation between the 
level densities around the bottom of the potential well, i.e., factor (i). For neutrons with A<- 
4 MeV, the value of |A| is even smaller: |A|<600keV. The largest deviations approaching 
2 MeV have been obtained for the neutron drip-line nuclei such as 122 Zr. Considering the 
analysis presented in this study, it has to be concluded that the excellent agreement found in 
Ref. |2D| (|A|~100keV) is fortuitous. According to our results in Table |, for nuclei discussed 
by these authors, i.e., 208 Pb and 208 114, the difference A is indeed very small. However, in 
other cases deviations are larger. 
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APPENDIX A: NEAR-THRESHOLD BEHAVIOR OF THE SEMICLASSICAL 

LEVEL DENSITY 



In the WK method, the semiclassical level density g sc (z) can be written as 

g sc (e) = #tf(» + 0-i (e), (Al) 

where ^tf(^) is the Thomas- Fermi (TF) level density and g~i(e) is the WK correction term 
(of the order of h^ 1 ). By employing Eqs. (|T3|) and (0), one can write g sc (s) as a derivative 
of the particle number with respect to e: 

= , ,-,<*) = -jA (A2) 



In Eq. (|A2|) , Atf is the TF particle number, 



4 /2M\* /•!•«(*) 



iVTF(e) = -^J y (e-y)^dr, (A3) 
while the TV 1 WK term is 

, T , n 1 /2mv r sc(£) 2j 

As usual, the classical turning point is defined by the relation V(r sc (e)) = e. 

The pathological behavior of g sc (e) close to the top of the potential barrier can be at- 
tributed to the singularity in the g~i(e) term. To examine this divergence, let us consider 
the integral 

1(e) = / -^=r 2 dr, (A5) 
•if y/e — V 

where it has been assumed that V(r) is an increasing function of r on an interval [f, r sc ], 
and f is a fixed radius (r < r sc ). 

By substituting x = V(r), 1(e) can be written as: 

= f -J0=dx, (A6) 
Jx y/e — X 



where x = V(r) and 

v (x)=r 2 ^— \ x = [r 2 — + 2r) \ x . (A7) 



,V 2 ^, / ,V" 



V \ V 

The singularity at x = e in the integrand in Eq. (|A6|) can be eliminated by performing a 
partial integration: 

1(e) = 2r](x)Ve - x + 2 / rf{x)y/e - xdx. (A8) 

Jx 



The first term in Eq. ( |A8|) does not cause any problems around the particle threshold and 
can be neglected in the following. Hence the behavior of g-i(s) around the top of the barrier 
is governed by the integral 

I'(e) « f 4^=dx. (A9) 
Jx \Je — x 
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1. Woods Saxon potential: neutron case 



For neutrons in a WS potential, the particle threshold appears at e=0. In the vicinity 
of the threshold, the potential energy can be approximated by 

V(r) « xexp(^-- — , (A10) 

and the classical radius is 

r(x) = f + a In f-^j . (All) 

In the limit e — > 0, the function 77 (x) in Eq. can be written as: 

r](x) = 2r(x) - r 2 (x)/a fa -r 2 (x)/a. (A12) 

Hence 

n'fo) ~ 2 ^ (A13) 

x 

and 

/■g ln(^) 

/'(e) « 2 / , w da;. (A14) 
xy£ — £ 

The above integral can be easily calculated. Around e=0, it behaves as 

In (I) /Jl, (A15) 
\x/ V x 

and this is the asymptotic behavior of g sc (e) around the neutron threshold. 

The behavior of g sc (e) for the neutrons in 120 Sn at is displayed in Fig. |8](a). It is 
seen that the semiclassical density diverges according to the law given by Eq. ( |A15| ). 

2. Finite potential barrier: proton case 

For potentials with finite barriers, such as the sum of WS and Coulomb potentials, the 
particle threshold appears at the top of the barrier, e—Vs- Around the barrier top, the 
potential energy can be expanded as 

V(r)f*V B -a^(r-r B ) 2 } (A16) 

where a = — V"{r B ). 

For e ~ Vb, the function 77 (x) in Eq. ( |A7| ) can be written as: 

r,(x) « (A17) 

re — r(xj 
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Consequently, 



and the leading term in I'(e) takes the form 



I'{e) -> / — rfs. (A19) 

(Vb — x) 6 l 2 y/e — x 



The above integral can be easily evaluated. Around e—Vb it behaves as 

1 



V B -s 



(A20) 



and this gives the asymptotic behavior of g sc (e) around the top of the Coulomb barrier. 

The behavior of g sc (e) for the protons in 120 Sn around e=Vb is displayed in Fig. §(b). 
It is seen that around the top of the barrier the semi-classical density diverges according to 
the law given by Eq. (A20). Of course, the WK contribution (A4) to the particle number 



also diverges when e — > Vb- This result is by no means surprising; the semi-classical approx- 
imation breaks down if the gradient of the potential at the turning point vanishes, and this 
happens precisely around the top of the barrier. 
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TABLES 



TABLE I. Shell correction 5E , the r.m.s. error a, and the Fermi level A calculated using the 
generalized Strutinsky method with continuum. The corresponding semi-classical quantities: shell 
correction 5E SC and Fermi level A sc are also shown together with the difference A=5E — 5E SC . All 
energies are in MeV. 









Neutrons 








Nucleus 


X TP 
OH/ 


a 


A 


X TP 

oE sc 


\ 

A sc 


A 


78 m: 

JNi 


-2.06 


0.183 


O C A 

-2.64 


-4.22 


-2.51 


i on 
1.39 




-7.iy 


0.100 


-9.63 


-6.82 


-9.77 


-0.37 


96 Zr 


0.24 


0.016 


-7.32 


0.82 


-7.37 


-0.58 


104 Zr 


6.57 


0.056 


-4.79 


6.48 


-4.71 


0.09 


106 Zr 


5.97 


0.039 


-4.23 


5.56 


-4.13 


0.41 


108 & 


5.76 


0.150 


-3.69 


4.94 


-3.57 


0.82 


no Zr 


4.49 


0.029 


-3.17 


3.45 


-3.05 


1.04 


122 Zr 


-4.61 


0.056 


-0.32 


-6.40 


-0.44 


1.79 


124 Zr 


-2.91 


0.052 


0.12 


-4.39 


-0.12 


1.47 


132 Sn 


-8.70 


0.023 


-4.50 


-8.94 


-4.42 


0.24 


146 Gd 


-10.09 


0.118 


-9.77 


-9.85 


-9.89 


-0.24 


208p b 


-11.37 


0.063 


-5.50 


-11.23 


-5.56 


-0.13 


114 


-8.44 


0.090 


-4.83 


-8.63 


-4.81 


0.19 








Protons 








Nucleus 


SE 


a 


A 


SE SC 


A S c 


A 


48 Ni 


-2.11 


0.084 


-0.16 


-1.94 


-0.08 


-0.17 


90 Zr 


1.96 


0.222 


-6.65 


1.45 


-6.84 


0.51 


100 Sn 


-7.31 


0.083 


0.72 


-7.01 


0.61 


-0.30 


132 Sn 


-6.02 


0.081 


-13.24 


-6.65 


-13.31 


0.63 


146 Gd 


5.27 


0.247 


-3.98 


4.51 


-4.12 


0.77 


180p b 


-7.72 


0.016 


-0.81 


-8.57 


-0.87 


0.85 


208p b 


-6.73 


0.028 


-7.16 


-7.33 


-7.22 


0.60 
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FIGURES 



FIG. 1. The distribution of Gamow energy eigenvalues, Wi in the (Re(u;), Im(w)) plane for 
(a) the stable nucleus 90 Zr (neutron eigenvalues), (b) neutron drip-line nucleus 122 Zr (neutron 
eigenvalues), and (c) proton-rich nucleus 180 Pb (proton eigenvalues). The contours L used in 
Eq. (|l7]) to calculate the smoothed level density in the generalized shell-correction method are also 
shown. Only the Gamow states with W{ lying above the contour are included in the leading term 
in Eq. (O). 



FIG. 2. The energy dependence of the "continuum particle number" (|l~8| ) for the neutrons in 
122 Zr (a) and the protons in 180 Pb (b) along the contours L of Fig. |l|(b,c). The energy dependence 
of N c along the path is very gradual. The fluctuations in N c can be attributed to the presence of 
near-lying Gamow states. 



FIG. 3. Shell correction for the neutrons in (a) 298 114, (b) 146 Gd, and (c) 90 Zr obtained using 
the generalized Strutinsky averaging procedure as a function of the smoothing range parameter 
7 for various orders of the curvature correction: p=8 (dotted line), p=12 (dot-dashed line), and 
p=16 (solid line). The continuum contribution to the level density has been calculated using the 
method described in Ref. 0. The gray line shows the result of the semi-classical Wigner-Kirkwood 
approach. 



FIG. 4. Comparison of the smoothed level densities calculated using the generalized Strutinsky 
method (SM) and the Wigner-Kirkwood method (WK) for (a) the neutrons in 146 Gd and (b) the 
protons in 208 Pb. The densities are normalized to the Strutinsky density g(e) calculated with the 
curvature correction p=10. Semi-classical level densities and Strutinsky level densities calculated 
with p=6 and 14 are shown by dotted, dot-dashed, and dashed lines, respectively. It is seen that 
the result of the Strutinsky smoothing is practically p-independent. The Fermi levels A (SM) and 
A sc (WK) are indicated, together with the value of the potential depth Vq. 



FIG. 5. Comparison of the smoothed level densities calculated using the generalized Strutinsky 
method (solid line, p=W variant) and the Wigner-Kirkwood method (dotted line) for the neutrons 
in 146 Gd (top) and the protons in 208 Pb (bottom). 



FIG. 6. Neutron shell corrections and Fermi energies as a function of N calculated in the SM 
and WK models using the single-particle potential of 90 Zr. The corresponding smoothed level 
densities are shown in Fig. 0. 



FIG. 7. Same as in Fig. || except for the neutrons in 90 Zr. The points at which the difference 
between SM and WK level densities, 5g, changes sign are marked by "A", "B", and "C". The 
oscillatory behavior of 5g is responsible for the oscillatory behavior of A as a function of particle 
number, as shown in Fig. |(| 
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FIG. 8. The divergent behavior of g sc (e) for the neutrons (a) and the protons (b) in 120 Sn 
around the particle threshold. It is seen that the semi-classical approximation breaks down in 
the vicinity of the threshold. The densities scaled according to Eqs. ( |A15[ ): 1000<7 sc (e)/[ln(e)/\/I] 
(where e = e/V Q ), and Eq. ( |A20|) : 1000c/ sc (V)/[i/e] (where e = (e-V B )/V ) are shown in the insets. 
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